function consJacBase ()
% Constraints and its "basic Jacobian" matrix, including Jacobian structure
%
global idx_gx idx_gy omgs TJ TJsum          % by initial.m
global Y3 Y3fun busCut rowCut colCut        % by Ymatrix.m
global uGalla Ia YxV                        % by Ymatrix.m
global prob nlp                             % nlp.func is modified
global rows cols JacBase locJx

calcBaseJacobian();                         % calculate JacBase

nlp.func.constraints = @constraintsFunc;
nlp.func.jacobianstructure = @() JacBase;

%% Functions
    function c = constraintsFunc (x, ~)                     % '~' for auto-differentiation
        % Constrains DAE
        nrow = nlp.nnode + 1; 
        %stateMat = reshape(x(1 : nrow*ncol), nrow, ncol);
        %odelhs = reshape(nlp.collocI.diffMat * stateMat, ...
        %    nlp.nnode *ncol, 1);                           % left hand side of ODE
        uG = x(nlp.uGidx); uL = x(nlp.uLidx);
        c = zeros(nlp.ncon, 1);
        odelhs = zeros(nlp.nnode * prob.nstate, 1);
        if isa(x, 'double') == 1          
            feval(Y3fun, uG, uL);                          % @ Ymatrix/updateY3 ()
        else                                               % !!! Auto-differentiation
            c = cons(c, x);                                % !!! Change class type
            odelhs = cons(odelhs, x);
            YxV = cons(zeros(2*prob.nB, 1), x); Ia = YxV;
            uGalla = cons(zeros(prob.nG, 1), x);
        end
        idxlhs = 1 : nlp.nnode; idxx = 1 : nrow;
        for ix = 1 : prob.nstate
            odelhs(idxlhs) = nlp.diffMat * x(idxx);
            idxlhs = idxlhs + nlp.nnode; idxx = idxx + nrow;
        end
        for inode = 1 : nlp.nnode
            dlt = x(nlp.dltidx(:, inode)); omg = x(nlp.omgidx(:, inode));
            V = x(nlp.Vidx(:, inode));     coi = x(nlp.coiidx(:, inode));
            dae = feval(prob.daefun, ...
                dlt, omg, V, coi, uG, uL);                 % @ modelOCP/Dae ()
            c(nlp.conidx(:, inode)) = -dae;                % -dae --> c
        end
        nlhs = nlp.nnode * prob.nstate;
        c(1 : nlhs) = odelhs + c(1 : nlhs);                % odelhs - oderhs
    end
    % Base Jacobian Matrix
    function calcBaseJacobian ()
        % auxiliary parameters
        nstate = prob.nstate; nG = prob.nG; nB = prob.nB;
        nnode = nlp.nnode;    onevec = ones(nnode, 1);
        [iD, jD, vD] = find(nlp.diffMat); nzD = length(vD);
        [iY, jY, vY] = find(Y3); nzY = length(vY);         
        % row numbers and column numbers
        rows = reshape(1 : nlp.ncon, nnode, prob.ndae);            
        nxvar = (nnode+1)*prob.nstate; nuvar = nnode*prob.nctrl; 
        xcols = reshape(1 : nxvar, nnode+1, prob.nstate);
        ucols = reshape(nxvar +1 : nxvar +nuvar, nnode, prob.nctrl);
        cols = [xcols(1:nnode, :) ucols];                         
        % determine Jacobian size
        nzJ = 2*nG * nzD + (9*nG + 1 + nzY + 2*prob.nparam) * nnode;
        nzJx = (5*nG + 4*prob.ncutG + 6*prob.ncutL) * nnode;
        iJ = zeros(nzJ, 1); jJ = zeros(nzJ, 1); vJ = zeros(nzJ, 1);  
        iJx = zeros(nzJx, 1); jJx = zeros(nzJx, 1);
        % Base Jacobian for: Ddelta - omgs(omega - 1) = 0
        %                    Domega - 1/TJ*(Pm - EXd*(Vgx*sin(delta)-Vgy*cos(delta))) = 0
        for istate = 1 : nstate                % diagonal of [states control; ODE]
            s = 1 + (istate -1) *nzD;    e = s - 1 + nzD;    idxJ = s : e;
            rowb = nnode * (istate - 1); colb = (nnode + 1)* (istate - 1);
            iJ(idxJ) = rowb + iD; jJ(idxJ) = colb + jD; vJ(idxJ) = vD;
        end                  % Non-zero<1>: 2*nG * nzD
        for idlt = 1 : nG                      % rows for Ddelta
            s = e +1; e = e +nnode; idxJ = s : e; vJ(idxJ) = -omgs; % d[F(Ddlt,...)]/d[omega]
            iJ(idxJ) = rows(:, idlt); jJ(idxJ) = cols(:, idlt + nG);             
        end                  % Non-zero<2>: nG * nnode
        for iomg = 1 : nG                      % rows for Domega: fake value=1
            r = rows(:, iomg + nG); rs = [r; r; r];
            cs = [cols(:, iomg); cols(:, 2*nG + idx_gx(iomg)); cols(:, 2*nG + idx_gy(iomg))];
            s = e +1; e = e + 3*nnode; idxJ = s : e;
            iJ(idxJ) = rs; jJ(idxJ) = cs; vJ(idxJ) = 1;             % d[F(Domg,...)]/d[delta, Vgx, Vgy]
            % Section <1> that will change [NZ = 3*nG * nnode]
            sx = 1+ (iomg - 1) * 3*nnode; ex = sx-1 + 3*nnode; idxJx = sx : ex;
            iJx(idxJx) = rs; jJx(idxJx) = cs;  
        end                  % Non-zero<3>: 3*nG * nnode
        % Base Jacobian for: 
        %                   Y * V - I = 0;
        iY = iY + 2*nG; jY = jY + 2 * nG;
        idxJ = e-nnode+1 : e;                                       % initialize idxJ
        for ie = 1 : nzY                       % Y matrix: uL and uG won't affect non-zero structure
            idxJ = idxJ + nnode; vJ(idxJ) = vY(ie);
            iJ(idxJ) = rows(:, iY(ie)); jJ(idxJ) = cols(:, jY(ie)); % d[F(V,...)]/d[V]
        end                  % Non-zero<4>: nzY * nnode
        e = idxJ(end);
        nzYx = 2 * prob.ncutG + 4 * prob.ncutL; 
        idxJx = ex-nnode+1 : ex;               
        for iyx = 1 : nzYx                     % Non-zero of changing part in Ymatrix(x)
            % Section <2> that will change [NZ = (2*ncutG + 4*ncutL) * nnode]
            idxJx = idxJx + nnode;
            iJx(idxJx) = rows(:, 2*nG + rowCut(iyx)); jJx(idxJx) = cols(:, 2*nG + colCut(iyx)); 
        end
        ex = idxJx(end);
        idxJ = e-2*nnode+1 : e; idxJx = ex-2*nnode+1 : ex;
        for ii = 1 : nG                        % I(delta): Igx = EXd *sin(dlt); Igy = EXd *cos(dlt);
            c = cols(:, ii); cs = [c; c];
            rs = [rows(:, 2*nG +idx_gx(ii)); rows(:, 2*nG +idx_gy(ii))];
            idxJ = idxJ + 2*nnode; vJ(idxJ) = 1;                    % fake value 1 for changing elements
            iJ(idxJ) = rs; jJ(idxJ) = cs;                           % d[F(V,...)]/d[delta]
            % Section <3> that will change [NZ = 2*nG * nnode]
            idxJx = idxJx + 2*nnode;
            iJx(idxJx) = rs; jJx(idxJx) = cs;
        end                  % Non-zero<5>: 2*nG * nnode
        for ip = 1 : prob.nparam                                    % d[F(V,...)]/d[uG, uL]
            cs = ones(2*nnode, 1) * (nxvar + nuvar + ip);
            rs = [rows(:, 2*nG +busCut(ip)); ...
                rows(:, 2*nG +nB +busCut(ip))];
            idxJ = idxJ + 2*nnode; vJ(idxJ) = 1;
            iJ(idxJ) = rs; jJ(idxJ) = cs;
            %  Section <4> that will change [NZ = 2*nparam * nnode]
            idxJx = idxJx + 2*nnode;
            iJx(idxJx) = rs; jJx(idxJx) = cs;
        end                  % Non-zero<6>: 2*nparam * nnode
        e = idxJ(end);
        % Base Jacobian for:  TJ' * dlt/TJsum - coi;
        %                     delta - coi;
        idxJ = e-nnode+1 : e; rowb = 2*nG +2*nB; rs = rows(:, rowb +1);
        for idlt = 1 : nG                                           % d[F(delta,coi)]/d(delta)
            cs = cols(:, idlt);
            idxJ = idxJ + nnode; vJ(idxJ) = TJ(idlt) * onevec;
            iJ(idxJ) = rs; jJ(idxJ) = cs;
            idxJ = idxJ + nnode; vJ(idxJ) = 1;
            iJ(idxJ) = rows(:, rowb +1 +idlt); jJ(idxJ) = cs;
        end                  % Non-zero<7>: 2*nG * nnode
        colb = rowb; cs = cols(:, colb+1);
        for icoi = 1 : 1+nG                                         % d[F(delta,coi)]/d(coi)
            idxJ = idxJ + nnode; 
            if icoi == 1, vJ(idxJ) = -TJsum; else vJ(idxJ) = -1; end
            iJ(idxJ) = rows(:, rowb + icoi); jJ(idxJ) = cs;
        end                  % Non-zeros<8>: (1 +nG) *nnode
        JacBase = sparse(iJ, jJ, vJ, nlp.ncon, nlp.nvar);
        locJx = sub2ind([nlp.ncon, nlp.nvar], iJx, jJx);
    end
end